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1 INTRODUCTION 

Several predictions obtained recently in the valence (quenched) approximation 
to the infinite volume, continuum limit of lattice QCD lie not far from experiment. 
For low lying hadron masses [1], valence approximation results are within 6% it 
8% of experiment. For decay constants [2] the valence approximation differs from 
experiment by increments ranging from 12% it 11% to 17% it 6%. Missing from these 
calculations, however, is an independent theoretical estimate of the error arising from 
the valence approximation. 

In the present article, we develop a systematic expansion for lattice QCD 
including the full effect of quark vacuum polarization. The leading term in this 
scheme is the valence approximation. If an infinite collection of higher terms is taken 
into account, full QCD is reproduced exactly. We then derive a formula which can be 
used to estimate the error in any vacuum expectation value obtained by truncation of 
this expansion to some finite number of terms. Our expansion assumes quarks occur 
in pairs with equal mass. 

In an exact treatment of QCD, virtual quark-antiquark pairs produced by a 
chromoelectric field reduce the field's intensity by a factor which depends both on 
the field's momentum and on its intensity. In the valence approximation this factor, 
analogous to a dielectric constant, is approximated by its zero-field-momentum zero- 
field-intensity limit [3]. The approximation which we consider here may be pictured 
as incorporating an inverse dielectric constant which is a sum of terms which pro- 
gressively more accurately reproduce the correct dependence of the inverse dielectric 
constant on field momentum and field intensity. 

In Section 2 we introduce definitions. In Section 3, we construct an expansion 
for the dependence of vacuum polarization on field momentum and field strength and 
derive an expression for the error in any vacuum expectation value arising from a 
truncation of the expansion. In Sections 4 and 5 two variations on the expansion 
of Section 3 are considered. In Sections 6, 7, and 8 we present an algorithm for 
evaluating the terms in these expansion. In Section 9, we describe a trial calculation 
using the expansion in Section 3. In Appendices A, B, C and D, we prove convergence 



of the sequence of vacuum expectation values arising from truncated forms of our 
expansions. Appendix E gives a calculation of a set of parameters needed by the 
algorithm in Section 6. 

The main motivation for the present work is to find a way of determining di- 
rectly from QCD the errors arising in valence approximation calculations of hadron 
masses and decay constants. This goal is accomplished in principle by the error for- 
mula for our expansion truncated after the leading term. A crucial question which 
we have not yet answered, however, is whether in practice the determination of va- 
lence approximation errors using our algorithm would be any faster than a direct 
comparison of valence approximation results with numbers found by the best present 
algorithms for full QCD. Whether or not the method we propose turns out to be 
useful in practice for quantitative error estimates, it appears to us that it does help 
provide a qualitative understanding of the mechanism underlying the relatively close 
agreement found in Refs. [1, 2] between valence approximation predictions and the 
real world. In particular, the results we present in Section 9 are clear evidence of 
at least one set of parameters and expectation values for which the main effect of 
vacuum polarization is absorbed by the dielectric constant implicit in the valence 
approximation. 

We are aware of two other strategies for evaluating the relation between the 
valence approximation and full QCD. The application of chiral perturbation theory 
to estimating the errors introduced by the valence approximation has been considered 
by several groups [4, 5, 6]. The limiting behavior at small quark mass of a variety 
of predictions of the valence approximation has been shown to be qualitatively dif- 
ferent from the behavior of full QCD. The quark mass below which these difficulties 
become quantitatively significant, however, appears to be well below the average of 
the up and down quark masses [7, 8]. For physical values of quark mass, several un- 
known parameters enter chiral perturbation theory predictions of the errors in most 
valence approximation results. Quantitative determination of these errors is there- 
fore not possible at present. Another method for evaluating of the effect of virtual 
quark-antiquark pair production on QCD predictions is discussed in Ref. [9]. This 
calculation uses a weak coupling expansion to leading order and is valid for small 



values of the gauge coupling constant and large values of the quark mass. The results 
we report in our trial calculation in Section 9 are qualitatively consistent with those 
described in Ref. [9]. Reviews of a variety of other recent valence approximation 
calculations are given in Ref. [10]. 

2 DEFINITIONS 

We consider Wilson's formulation of Euclidean QCD on some finite lattice. A 
lattice gauge field consists of an assignment of an element u[xi^ X2) of the fundamental 
representation of SU{3) to each oriented nearest neighbor pair of sites (a;i,a;2) with 
the usual restriction that u{xi,X2) is the adjoint u{x2,xiy . 

Define the Hilbert space JF to consist of complex valued functions / of the 
lattice gauge fields with finite value of the norm 

ll/ir = C-\j d^\f\' exp{S), (2.1) 

C = I djjL exp{S). 

The inner product on T is 

(/,/') = C\f dfif* f exp{S), (2.2) 

( = I djjL exp{S). 

Here S is some real-valued function of the field which is bounded in absolute value. 
Several different possible choices of S will be discussed in Sections 3 and 4. A linearly 
independent basis for T consists of the collection {/j} of all possible products of 
matrix elements of irreduceable representations of SU{'i) including exactly one matrix 
element for each link, with links differing only by a flip of orientation identifled. 
Distinct fi are then orthogonal with respect to the inner product of Eq. (2.2) for S 
of 0. Eor convenience in Appendix C, we choose the fi to be normalized with respect 
to the inner product with S of 0. 

Let di be the sum over all links of the dimension of the SU{'i) representation 
assigned to that link by fi. We assume the sequence {fi} is ordered in such a way 



that di is a nondecreasing function of i. Applying a Gram-Schmidt process to {fi} 
using the inner product of Eq. (2.2) for some nonzero choice of S gives an orthonormal 
basis {fi} for JF. 

Although the expansion to be constructed in Section 3 can be defined using 
only JF, for purposes of constructing an algorithm to evaluate this expansion it is 
slightly more convenient to work with the subspace "H of JF which is invariant under all 
lattice translations, rotations, refiections, gauge transforms and complex conjugation. 
Let hi be the projection of fi onto Ti. Since rotation, translation, refiection, gauge 
transformation and complex conjugation leave the value of di unchanged, hi will be 
a linear combination of a collection of fj all of which have dj equal to di. Most hi 
obtained in this way will be linearly dependent on the set of hj with j < i. Working 
upwards from z of 0, we eliminate any hi which is dependent on surviving hj with 
j < i. A Gram-Schmidt process on the surviving sequence gives an orthonormal basis 

{k} for n. 

Typical vectors in Ti are the function with value 1 for all gauge fields and the 
Wilson plaquette action 

P = 22 tr[u{xi,X2)u{x2,X3)u{x3,X4)u{x4,Xi)], (2-3) 

(xi,...X4) 

where the sum is taken over all oriented plaquettes (xi, . . . X4) consisting of sequences 
of four successive nearest neighbors, with sequences related by a cyclic permutation 
identified. Any sum of traces of products of u[x^y) over all rotations, translations, 
refiections and order reversals of some closed path gives yet another element of Ti. 
The basis vector ho is the function with constant value 1, and hi is the normalized 
projection of P orthogonal to ho. The basis vectors /i2, /^s, and /i4, are each found by 
continuing the Gram-Schmidt process on the three different sums of traces of products 
of u[x^ y) along one of the three distinct shapes of closed paths consisting of six lattice 
links. 

We now define the lattice vacuum expectation value. We assume quarks occur 
in degenerate pairs for some set of masses strictly greater than 0. Let M be Wilson's 
coupling matrix among half the quark fields, one from each degenerate pair. We 



impose periodic boundary conditions. For any function of the gauge fields G with 
bounded absolute value, a regulated form of the vacuum expectation value obtained 
after integrating out quark fields is 

<G>R = Z-^ f clfiGdetiM^M + R)exp[-P], 
J 6 

/5 



Zr = I diidet{M^M + K)exp{-P), (2.4) 



where (3 is 6/(/^ for bare gauge coupling constant (/, /i is the product of one copy of 
SU{3) Haar measure for each link variable on the lattice, and R a small nonnegative 
parameter. The extension of Eq. (2.4) to vacuum expectations of products of quark 
and antiquark fields is not needed for the present discussion and will be omitted for 
simplicity. Since the Wilson coupling matrix M obeys 

det{M) = det{M^) (2.5) 

the expectation < G >o is the usual vacuum expectation of lattice QCD. For any 
bounded G, < G >r approaches < G >o as i? goes to 0. 

We introduce the regulator R in the definition of < G >r to provide a mathe- 
matically convenient rule for handling rare gauge configuration on which M becomes 
singular in the valence approximation. Monte Carlo valence approximation calcula- 
tions often find averages of quantities involving M~^ at values of the quark mass for 
which some configurations exist, such as all link variables close to the identity matrix, 
for which M has eigenvalues arbitrarily close to 0. These configurations are not en- 
countered in practice because their total weight within the path integral is extremely 
small. It is generally believed that for any positive choice of quark masses, the total 
valence approximation measure of configurations with minimal M^ M eigenvalue be- 
low 0{m^) goes to zero very rapidly in the limit of large lattice volume with lattice 
spacing held fixed. The expansion to be considered below will be done with some 
nonzero value of R much smaller than 0{m^). After taking a limit of infinite volume 
of any vacuum expectation, a limit of zero R should leave the result essentially un- 
changed. For notational simplicity, the R subscript will be deleted from < G >r and 
Zr in the following. 



The regulation parameter R is needed only for the Wilson quark coupling 
martix. For Kogut-Susskind quarks the spectrum of M^ M is bounded from below by 
m'^ for all gauge configurations. 

3 EXPANSION 

For any choice of iS, bounded in absolute value, in Fq. (2.2), the space Ti 
can be used to construct an expansion for logdet{M^ M + K). Since M is a finite 
matrix with matrix elements bounded uniformly over all gauge fields, the spectrum 
of M^ M + i? is bounded from above by some constant A. Thus det{M^M + K) has a 
finite value of the norm defined by Fq. (2.1) and is in JF. In addition det{M^ M + K) 
is real- valued and rotation, translation, refiection and gauge invariant. It is therefore 
inn. 

Since the spectrum of M^ M + i? is bounded from below by R, det{M^ M + K) 
is bounded from both above and below. Thus || log o?et(M^M + K)\\ is finite and 
logdet{M^ M + K) is also in Ti. Using the orthonormal basis {hi} of Ti, we therefore 
have the convergent expansion 

logdet{M^M + R) = ^a,A„ (3.1) 

i 

a, = {k,logdet{M^M + Rj). (3.2) 

An algorithm for the numerical evaluation of the coefficients a^ is presented in Sec- 
tion 6. 

For any G bounded in absolute value, an approximation to < G > can be 
obtained by combining Fqs. (3.1) and (3.2) with Fq. (2.4). The expectation value 
defined by Fq. (2.4) can be reexpressed 

<G> = Z-^ J dfiGeXplLn + Qn + ^P], 

Z = dijexp[Ln + Qn + -P], (3.3) 



where n is some positive integer and the partial sum L^ and remainder Q^ are 

n 

Ln = ^a,h,, (3.4) 

8 = 

Qn = logdet{M^M + R)-Ln, 

= ^a,h,. (3.5) 

On 

As n becomes large, Eq. (3.1) implies Qn approaches in T. Thus it appears reason- 
able to try to approximate Eqs. (3.3) by omitting Qn- We obtain 

< G >n = Z^^ diiGexp[Ln + -P], 

Zn = dfiexplLn + -P]. (3.6) 

The expectation < G >o is pure QCD with the quark determinant simply removed 
and no shift in (3. The expectation < G >i is the valence approximation including a 
shift in (3. For any G bounded in absolute value, the approximate expectation < G >„ 
approaches < G > as n becomes large. A proof is given in Appendices A, B, C and 
D. 

The asymptotic expansion to leading order in Qn for the error in < G >„ is 
the correlation 

<G> - <G>n = <{Qn- <Qn>n){G- <G>n)>n- (3.7) 

Combining Eq. (3.7) and the Cauchy-Schwarz inequality we can obtain the bound 

{<G>-<G >nY < < [Qn- < Qn >nY >n< (G- < G >nY >n • (3.8) 

If the estimate for < G > — < G >„ in Eq. (3.7) turns out to be small in 
comparision to < G >„, the approximation < G >„ and the error estimate < G > 
— < G >n should both be reliable. Conversely if the error given by Eq. (3.7) is 
significant in comparision to < G >„, neither the approximate expectation < G >„ 
nor < G > — < G >n can be considered reliable. Seen through the eyes of weak 
coupling perturbation theory, Eq. (3.7) gives the one quark loop approximation to 
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the error < G > — < G >„. Thus the rehabihty of this error estimate is not directly 
related to the rate of convergence of the sequence of < G >„ as n is made large. 
Higher order terms in the loop expansion for < G > — < G >„ will be discussed in 
slightly more detail in Section 5.1. 

To try to minimize the valence approximation error <G> — <G>i, we now 
take S in Eq. (2.2) to be the effective action entering Eq. (3.6) for < G >i, 

S = L, + ^P, 

D 

= « + ^^, (3.9) 

with a and /?' given by 

a = ao -\- 



ai < P >i 



'< P- < P>i]2 > 

/?' = /3 + A/3, 
A/3 = ^""^ ^. (3.10) 

'<P-<P >i]2 > 



If S is held fixed in the inner product (...,...) and in the definition of the vacuum 
expectation < . . . >i but Oq and ai are varied in Qi, the values of Oq and ai given by 
Eq. (3.2) using (...,...) with S of Eq. (3.9) minimize < {Qi— < Qi >iY >i- Thus 
by Eq. (3.8) the error {< G > — < G >i)^ should tend to be minimized by the choice 
Eq. (3.9). 

Eor any choice of /3, the valence approximation value fi' is to be found by 
solving 

/3' = /3 + A/3(/3'). (3.11) 

With any reasonable number of fiavors of quarks, less than 10 for example, it is eas- 
ily confirmed numerically that A/3(/3') is a monotone increasing function of fi' with 
derivative significantly less than 1. Eor two fiavors of quarks with k of 0.1600 con- 
sidered in Section 9, we found that A/3 rises from below 0.1 at /?' of to below 0.3 
at fi' near 6.0. Eq. (3.11) can then be solved by a fixed point iteration taking fi as 
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an initial value of /?'. Alternatively, if the valence approximation /?' is chosen first, 
Eq. (3.11) gives directly the corresponding (3 for full QCD. This procedure will be 
adopted in the example to be discussed in Section 9. 

It may be useful to mention here that the hopping constant expansion for 
logdet{M^ M) expresses this quantity as a linear combination of Wilson loops for- 
mally similar to Eq. (3.1) and can be used to obtain an approximation to < G > 
similar to Eq. (3.6). In two crucial ways, however, the expansion of Eq. (3.1) dif- 
fers from the hopping constant expansion, and approximation Eq. (3.6) differs from 
the corresponding approximation using the hopping constant expansion. Eirst, the 
validity of expansion Eq. (3.1) and the accuracy of approximation Eq. (3.6) are not 
restricted to the range of large quark mass to which the hopping constant expansion 
and its related approximation apply. As we have already shown and as will be dis- 
cussed in Appendices A, B, C and D, expansion Eq. (3.1) and approximation Eq. (3.6) 
apply as long as the quark mass is greater than 0. Second, even for values of quark 
mass at which the hopping constant expansion does converge, Eqs. (3.1) and (3.2) 
differ from the hopping constant expansion by an infinite rearrangement. That is 
each term which appears in Eqs. (3.1) and (3.2) is a linear combination of an infinite 
set of the terms appearing in the hopping constant expansion, and vice versa. 

A wide range of other possible expansions and approximate vacuum expec- 
tation values similar to Eqs. (3.1) and (3.6), respectively, can be constructed by 
changing the choice of S and expansion basis {hi}. Some of these will be discussed in 
Section 4. Yet another class of possibilities, which we will not discuss here in detail, 
is to choose the a^ in Eq. (3.1) to force the first order shift in Eq. (3.7) to zero for a 
particular G or set of G, such as the pion propagator or the rho propagator. 

4 OPTIMIZED EXPANSION 

In Section 3, we showed that the leading contribution to a valence approxima- 
tion error {< G > — < G >iy will tend to be minimized by choosing the effective 
action S in Eq. (2.2) to be the effective action Li -\- ^P occurring in the valence ap- 
proximation vacuum expectation < ... >i. Since the effective action contribution Li 
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depends on the choice of iS, we are then led to a nonhnear equation Eq. (3.11) to be 
solved to find S. A similar recursive definition of iS as iv„ + ^P will tend to minimize 
the leading contribution to the error (< G > — < G >n)^- This choice gives a more 
complicated version of Eq. (3.11) to be solved for the coefficients entering S. 

In Section 3 we also argued, however, that the function A/3(/3) defined by 
Eq. (3.10) has a small derivative. This observation can roughly be rephrased as 
saying that in the neighborhood of our choice of iS, the valence approximation effective 
action coefficients vary only slowly with respect to changes in the S entering Eq. (2.2) 
defining (...,...). We would expect similarly weak dependence on S of the effective 
action L^ + ^P of higher order vacuum expectations < . . . >„. Thus the improvement 
obtained in any (< G > — < G >n)^ by using a complicated S in Eq. (2.2) in 
place of the S of Eq. (3.9) may not be very large. Similarly, an S which is a close 
approximation to P„ + ^P may provide very nearly all of the decrease available in a 
typical {< G > — < G >n)^ by optimal choice of S. Rather than trying to solve a full 
recursive equation to optimize S for < . . . >„, a choice close to optimal is probably 

Sn = L^ + §P, (4.1) 

D 

with Ln itself taken from Eq. (3.4) using S of Eq. (3.9). A choice yet closer to optimal 
would be to apply Eq. (4.1) iteratively taking P„ from an inner product defined using 

5 YET ANOTHER EXPANSION 

If < G > in Eq. (3.3) is expanded as a power series in Q^ < G >„ of Eq. (3.6) 
gives the constant term, and < G > — < G >„ of Eq. (3.7) gives the linear term. 
Continuing to expand < G > in powers of Qn yields an infinite series. It is likely that 
the resulting series can not be summed and is only an asymptotic expansion in small 
Qn- An alternate expansion for < G > as a ratio of power series in Q^ however, 
can be proved convergent. The leading approximation to this ratio also reduces to 
< G >n plus the first correction in Eq. (3.7). 
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Expanding the two exponentials which appear in Eq. (3.3) as power series in 
gives 

Y 



<G> 






m\ 



The series in Eq. (5.1) are both convergent since Qn, for any choice of n, is bounded 
in absolute value by a constant independent of the gauge field. This bound on Qn 
holds since both log o?et(M^M + i?) and Ln in Eq. (3.5) defining Qn obey such bounds. 
The algorithm to be described in Section 6 can be adapted to evaluate each of the 
individual terms in the series in Eq. (5.1). Evaluating terms beyond ?7i of 1 using this 
method, however, is probably very time consuming. 

The terms of order m in Eq. (5.1) give, roughly, the m quark loop contribution 
to the error in < G >„. The truncation of the sums in Eq. (5.1) yields a two parameter 
family of possible approximations to < G >. Eor n and m both 1, < Qi >i vanishes 
so that Eq. (5.1) is equivalent to adding to < G >i the error estimate of Eq. (3.7). 

6 TRACE LOG ALGORITHM 

The quantity log o?et(M^M + R) needed for Eq. (3.2) obeys the identity 

log det{M^M + R) = tr log{M^M + R). (6.1) 

We now consider an algorithm for finding tr log(M^M + K). The algorithm exploits 
properties of the Chebyshev polynomials. Combined with Eq. (3.2), this algorithm 
gives the coefficients a^. As discussed in Section 2, we will assume the quark masses, 
lattice volume, and Monte Carlo ensemble size have been chosen in such a way that 
for all gauge configurations encountered the minimal eigenvalue B of M^ M lies well 
above R. The effect of R in trlog(M^M + K) is then negligable, and we omit R in 
the remainder of the paper. 
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To evaluate tr log (M^M) we begin by generating an ensemble of gaussian ran- 
dom complex-valued pseudo-quark fields (j)i{x)^ where z is a multi-index ranging over 
all combinations of quark spin, color and flavor and x ranges over lattice sites. For 
each i and x we choose (j)i{x) to be an independent random variable such that the 
average over this ensemble << ... >> gives 

« (f),{x) cf)j{y) »=0, 
<< (j)*{x) (j)j{y) »= % 5,y. (6.2) 

We then have 

trlog{M^M) = <<((</), log(MtM)</))) >>, (6.3) 

where ((...,...)) is the inner product on the vector space of pseudo-quark fields 

ix 

Finding the inner product of two such vectors requires a comparatively small amount 
of arithmetic. The problem of evaluating the trace trlog(M^M) is thus reduced to 
finding log{ M^ M)(f) for a large ensemble of (/>. 

For the evaluation of log{ M^ M)(f) we combine properties of the Chebyshev 
polynomials with the restriction that the eigenvalues of M^ M lie between upper and 
lower bounds A and 5, respectively. Define the operator Y and the parameter e to 
be 

Y = — ^, (6.5) 

e = ^. (6.6) 

In Appendix F, we will show that for any n greater than 1 there are a set of coefficients 
6i, such that for any number y between e and 1 

logy = ^6,7;*(^) + <^logy, (6.7) 

8 = 

\8\ < 2exp{-2ny^), (6.8) 
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where the T* are Chebyshev polynomials. For large values of n, the inequality of 
Eq. (6.8) is nearly saturated. Since K is a self-adjoint operator with all eigenvalues 
between e and 1, for any vector (j) 

n ] —Y 

log(y)</> = Y.b,T*{-—)cf^ + Slog{Y)cf^, (6.9) 

8 = 

with 5 bounded according to Eq. (6.8). 

An iterative algorithm to evaluate the sum in Eq. (6.9) can be obtained from 
the recursion relation 

T*^,{z) = {4z-2)T*{z)-T*_,{z) (6.10) 

and initial values 

To{z) = 1, 
Tl^{z) = 2z-l. (6.11) 

Define the sequences Ci and Di for < i < n hy 

a = T*iY)^ + c,i—j— a_i - A-i)</', 

A = T*_,{Y)(j) + c,a.,(j), (6.12) 

with initial values 

Co = T*{Y)<P, 

Do = T*_,{Y)(j). (6.13) 

The coefficients q in Eq. (6.12) are found from the bi in Eqs. (6.7) and (6.9) by 

Q = ^^. (6.14) 



-'n — t 



Eqs. (6.10 - 6.14) imply that boCn gives the sum in Eq. (6.9) and is therefore an 
approximation to log[Y)(j) with relative error less than \5\. 
The final result, by Eqs. (6.3) and (6.5), is 

trlog{M^M) = «((</), 6oC„))»+(/Mlog(A), (6.15) 

where o?m is the dimension of the matrix M. 
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7 MAXIMUM AND MINIMUM EIGENVALUES 

Since for small values of 8, the number of iterations required to obtain a fixed 



value of 5 in Eq. (6.7) becomes a linear function of J A/ B^ the optimal choices for A 
and B become Xmax and Xmim respectively. 

An efficient algorithm [12] for estimating the maximum and minimum eigen- 
values of M^ M uses the Lanczos method to construct a tridiagonal approximation 
to M^M. Define the sequences of real numbers «!,...«„ and /3o,.../3m and the 
sequences of pseudo-quark fields go, • • • ?m and Tq, . . . rm by 



87 



q^+i = {(3,) Vj 

a,+i = ((g,+i,MtMg,+i)), 



r 



A+i = ^((r,+i,r,+i)), (7.i; 



with /3o of 1, go identically 0, and Tq a randomly chosen pseudo-quark field with norm 
1. Here m is the number of distinct eigenvalues of M^M. 

It can be shown that the sequence of pseudo-quark fields gi , . . . g^, generated 
by Eq. (7.1) is orthonormal, and the space spanned by these vectors is invariant under 
the action of M^ M. The space spanned by gi, . . . g^, is smaller than the whole space 
of pseudo-quark fields only if one or more of the eigenvalues of M^ M is degenerate. 
In the basis gi, . . . g^,, M^ M is tridiagonal. All matrix elements of Tij 

T,, = {{q,,M^Mq,)) (7.2) 

vanish except 



T.-u -- 


= A-i, 


r„- = 


= a», 


Ti+ii - 


= A- 



(7.3) 

Each distinct eigenvalue of M^ M occurs exactly once as an eigenvalue of the 
matrix of Eq. (7.3). As n grows, the maximum and minimum eigenvalues AJ^^^ and 
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^min of the submatrix T^^ with 1 < «, j < n approach the true maximum and minimum 
eigenvalues Xmax and Xmin of M^ M with errors faUing exponentially in n. 

To extract AJ^^^ and X^i^ from T", define the polynomial p"(A) to be the 
determinant 

p"(A) = c/et(r" - A/"), (7.4) 

where /" is the n X n identity matrix. The p"(A) can be calculated from the iteration 

/(A) = 1, 
p"+i(A) = K+i-A)p"(A)-(/3„_i)V"nA). (7.5) 

The eigenvalues T" are the zeros of p"(A). Thus we wish to find the largest 
and smallest of these zeros. It can be shown [12] that the number of zeros of p"(A) 
lying below some A is given by the number of sign changes in the sequence 

p'{X),p'{X),...p-{X). (7.6) 

This relation can then be used to guide a search for the maximum and minimum 
zeros of p"(A). 

8 PRECONDITIONER 

With A and B given by X^ax and Xmim respectively, the amount of arithmetic 



required by the algorithm of Sect. 6 is proportional to y X^ax/ Xmin- We now show how 
the calculation of trlog(M^M) can be converted into the calculation of trlog(A^^A^) 



for an operator N with a smaller value of yX^ax/ Xmin- This change also tends to 
decrease the number of pseudo-quark fields needed for a reliable evaluation of the 
trace. 

The expression tr log (M^M) is gauge invariant since it is log det{M^ M) and 
det{M^M) is gauge invariant. Prior to evaluating trlog(M^M) we can therefore 
transform to a lattice transverse gauge, defined to give a local maximum of the sum 
over all nearest neighbor x and y 

J2tr[u{x,y)]. (8.1) 



xy 
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Using, for example, the algorithm described in Ref. [1], the number of arithmetic 
operations required for gauge fixing is relatively small in comparison to the arithmetic 
needed to find log{ M^ M)(f) for an ensemble of random pseudo-quark fields (/>. 

Now define Mq to be the fermion coupling matrix with hopping constant ko 
and all u[x^y) equal to 1. Since M has been transformed to a smooth gauge, if the 
bare gauge coupling constant g is made small and ko is chosen optimally, we expect 
Mo to be an approximation to M. Thus the preconditioned operator N 

N = (Mo)-'M, (8.2) 

should be closer to the identity than is M. In particular, J \max I ^min for N'^ N should 
be smaller than it is for M^ M. Using the preconditioned operator we have the relation 

trlog{M^M) = trlog{N^N)+trlog{M^oMo). (8.3) 

The additional term trlog(MjMo) required to find trlog(M^M) by Eq. (8.3) 
does not depend on the gauge configuration and needs to be calculated only once. On 
the other hand, the operator Mq is diagonal in momentum space. Thus fast fourier 
transforms provide an efficient way to carry out the multiplication by negative powers 
of MqMq needed to determine log{ N^ N)(f) by the algorithm of Sect. 6. A rough guess 
might be that \/ Xmax / ^min for N^ N will go to a constant if g is made small, while 



^max / ^min for M^ M wiU progrcssivcly grow. Thus it seems plausible that for small 
enough g the additional cost of fourier transforms required to apply the algorithm of 
Sect. 6 to the preconditioned operator will be more than made up for by the decrease 
in y^max/^min ^ud Corresponding decrease in the number of iterations of Eq. (6.12). 
At least for the set of parameters at which we run the algorithm in the example 
describe in the next section, this expectation turns out to be correct. 

9 EXAMPLE 

As a first test, we applied the algorithms of Sects. 6 - 8 to QCD with two 
fiavors of quarks both with k of 0.1600 on a 6'^ lattice. We calculated a variety of 
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vacuum expectation values using Eq. (3.6) with n of 0, with n of 1, which is the 
valence approximation, and then found the corresponding error estimates for n of 1 
from Eq. (3.7) and corrected vacuum expectations including these error estimates, 
equivalent to Eq. (5.1) with n of 1 and m of 1. These results were all compared with 
direct calculations in full QCD. 

All calculations in this section were done on an IBM RS/6000 workstation 
sustaining approximately 10 Mflops. The final production runs with our algorithm 
required about one week of machine time. The final comparison calculation with full 
QCD took about another one week. Another month or so of machine time was spent 
checking and exploring. 

As discussed in Section 3, rather than fixing /3, we fixed /?' given by the sum 
(3 + A/3. Our first task was then to calculate A/3. Erom /?' and A/3 we determined (3 
for full QCD. Eor /?' we chose the value 5.700. Then < . . . >i becomes simply a pure 
gauge vacuum expectation with pure gauge /?' of 5.700, and < . . . >o is a pure gauge 
vacuum expectation with the same (3 as used in full QCD, 5.700 — A/3. 

To evalulate the expectations < . . . >o and < ... >i, ensembles of pure gauge 
configurations were generated using the Cabbibo-Marinari-Okawa algorithm. Eor 
< . . . >o /? is comparatively small and we were not concerned with obtaining great 
precision. We found it sufficient to use 100 configurations with 100 sweeps to equili- 
brate and 100 sweeps between successive pairs. Eor < ... >i, however, we used 160 
configurations with 1000 sweeps to equilibrate and 1000 sweeps between successive 
pairs. Eor all of the quantities for which < ... >i was measured, we found 1000 
sweeps to be more than sufficient to produce equilibrium values and to decorrelate 
successive values. The full QCD results were found using the hybrid Monte Carlo 
algorithm. Hamiltonian trajectories were generated using the algorithm of Ref. [11], 
which is faster than leap-frog by about a factor of 2. Vacuum expectations were 
taken over an ensemble of 250 accepted trajectories each of length 0.25 time units, 
with 150 trajectories at (3 of 5.44 followed by 25 trajectories at (3 of 5.439 used to 
obtain equilibrium. 

To tune the algorithm of Sect. 7 for < ... >i, we began by evaluating on 
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each gauge configuration tlie ratio J \max I ^min for M^M. For tlie preconditioned 
operator N^N, we evaluated tliis ratio for a range of ko and found tlie ko which 
minimizes \/\max/ ^min- Results are shown in Table 1. The total work required to 
calculate log{M^ M)(f) is expected to be about 35% greater than the work required to 
find log{ N^ N)(f) to the same accuracy. Trial calculations were consistent with this es- 
timate. In the remainder of this section we therefore consider only the preconditioned 
operator and Eq. (8.3) to find trlog(M^M). 

Using the preconditioned operator N with the optimal A;o, we then calculated 
the expectation values 

El = <trlog{N^N)>i 

E2 = < [trlog{N^N)- < trlog{N^N) >i][P- < P >i] >i, (9.1) 

for a range of different choices of the number of iterations of the Chebyshev algorithm, 
Eq. (6.12). The results are shown in Table 2. The averages in Table 2 were found using 
a collection of 16 gauge configurations and 20 random (j) for each configuration. Eor 
each configuration, we first calculated \/\max/ ^min- The number of iterations of the 



Chebyshev algorithm, Eq. (6.12), was then chosen to be proportional to y Xmax / ^min ■ 
The value of rich shown in Table 2 is the number of iterations which would result for 



the average over all 160 configurations < yXmax/^min >i- As rich is increased above 
50, the change in the two measured expectations shown in Table 2 is significantly less 
than the statistical errors in the expectations we found with our final, full ensemble 
of gauge configurations and (j). Eor the remaining calculations, we chose rich to be 50. 
Using the final ensemble of 160 gauge configurations and a range of values of 
the number n^ of random (j) for each gauge configuration we calculated Ei and E2 
of Eq. (9.1). We also evaluated the dispersion < [P— < P >i]^ >i. We obtained 
the value (5.68 ± 0.61) x 10^ Erom E2 and < [P- < P >i]^ >i we then found A/3. 
The results are shown in Table 3. As expected, the values we obtained are consistent 
within errors as n^ is varied while the size of the errors themselves tends to fall as n^ 
increases. The optimal choice of n^ producing the smallest statistical uncertainty in 
A/3 for a fixed amount of computation can be shown to be roughly 100. 
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The error bars on the numbers in these tables are statistical, found by the 
bootstrap method. From the ensemble of 160 data sets, each consisting of a gauge 
configuration and an associated collection of n^ random (/>, we randomly chose 160 
new data sets to generate a bootstrap ensemble. On each bootstrap ensemble we 
then found Ei. In this way 100 bootstrap ensembles and 100 values of Ei were found. 
From these we evaluated the difference between the value of Ei larger than all but 15 
results, and the value of Ei smaller than all but 15 results. Half of this difference is 
shown as the statistical error. The errors for E2 and A/3 were found similarly. In the 
calculation of the error for A/3, < [P— < P >iY >i was calculated independently on 
each bootstrap ensemble and used to determine the corresponding bootstrap value 
for A/3. 

The most reliable value for A/3 is 0.261 it 0.014, obtained with n^ of 140. 
Our algorithm then predicts that expectation values in full QCD with two flavors of 
quarks, k of 0.1600 and /3 of 5.439 will agree with < . . . >i at /3' given by /3 + A/3 
which is 5.700 ± 0.014. 

Figure 1 shows vacuum expectations of a collection of different Wilson loops. 
All loops are rectangular with dimensions as shown except for 61 consisting of the 6 
link boundary of a pair of orthogonal plaquettes joined on an edge, and 62 consisting 
of the 6 link boundary of three orthogonal plaquettes joined to form half the surface 
of a cube. The normalization of each loop is 

-tr[u{xi, X2)U{X2, X3) ... u{xn, ^i)], (9.2) 

so that if all link matrices u[x^ y) were given by the identity matrix, all loops would 
become 1. Since the error bars on all points in Figure 1 are smaller than the symbols, 
in Table 4 we also give numerical values. Fach loop expectation shown in the flgure 
is obtained in four different ways. Boxes indicate Fq. (3.6) with n of and /3 of 5.439. 
Triangles represent the valence approximation, Fq. (3.6) with n of 1 and /3 + A/3 of 
5.700. Circles give the valence approximation, Fq. (3.6) with n of 1 and /3 + A/3 of 
5.700 but corrected by the error estimate of Fq. (3.7), equivalent to Fq. (5.1) with n 
of 1 and m o{ 1. Plusses show full QCD with /3 of 5.439 and two flavors of quarks 
with k of 0.1600. For all but the 1x1 loop, clear differences can be seen in Figure 1 
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between data points with n of and (3 of 5.439 and points for the full theory with two 
flavors of quarks and (3 of 5.439. For the largest loops these two results differ by as 
much as a factor of f 0. For most loops, nearly all of this shift is correctly reproduced 
by the valence approximation, n of 1 and (3 + A/3 of 5.700. For the largest loops 
the valence approximation reproduces most of the shift due to vacuum polarization, 
but falls noticeably below the results of full QCD. When the leading error estimate 
of Fq. (3.7) is added to the valence approximation, however, in all cases agreement 
with full QCD is found within statistical errors. Thus Fq. (3.7) appears to be quite 
reliable as a error estimate for the valence approximation for the data considered in 
Figure 1. 

A comparison between the numerical efficiency of our algorithm and that of 
the hybrid Monte Carlo method shows that our method is more efficient for one 
goal but less efficient for another. For the hybrid Monte Carlo method, we grouped 
successive configurations into bins and found the averages over each bin. We then 
evaluated the dispersion in the full ensemble averages by applying the bootstrap 
method to the binned ensemble. As the size of the bins used in this calculation is 
made larger, successive bin averages become statistically independent so that the 
dispersion predicted in this way for the full ensemble average becomes independent 
of bin size. We determined that about 16 hybrid Monte Carlo trajectories are needed 
to produce a new statistically independent configuration. On the other hand, since 
our method uses only a pure gauge updating algorithm, it is relatively inexpensive to 
guarantee the statistical independence of each successive configuration on which an 
ensemble of random (j) is constructed. The number of arithmetic operations required 
to generate one independent configuration by hybrid Monte Carlo turns out to be 
sufficient to generate about 7 new configurations and (j) ensembles by our method, if 
the optimal value of n^ is chosen. 

To obtain a first estimate of a vacuum expectation by either method we might 
require at least 16 statistically independent configurations be used. Using fewer than 
perhaps 16 configurations, it is difficult to determine with much confidence whether 
any independent, equilibrium configurations have been generated. According to this 
assumption, a first estimate of a vacuum expectation by hybrid Monte Carlo requires 
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about 7 times as much work as by our method. The statistical uncertainty in the 
hybrid Monte Carlo result found in this way, however, is significantly smaller than that 
determined by our method. For the data in Figure 1 and Table 4, with approximately 
equal time spent on full QCD and on the first correction to < . . . >i expectations, 
the statistical errors on the full QCD data are about a factor of 3 smaller than the 
statistical errors on the sum of < ... >i and its first correction. What the relative 
performance of the two different algorithms would be for values of lattice spacing, 
lattice volume and quark mass closer to the continuum limit, we do not yet know. 
As a further check of our method, we have calculated the expectation values 

^3 = < {tr[log{N^ N)]- < tr[log{N^N)] >i}2 >i, 
E, = <{Qif>i, 

= < [log{N^N) - L,]' >i, (9.3) 

= < {trlog(N^N)- < trlog(N^N) >i -^[P- < P >i]}' >i, 

6 

for Qi of Fq. (3.5), Li given by Fq. (3.4) and A/3 of 0.261. A small value of £"4 in 
comparison to E^ suggests that Li is a good approximation to tr log(A^^A^). 
To find these expectations, we evaluated 

E, = —<[J2iilogiN^N)^,,logiN^N)^i))]>^, (9.4) 

^^ .=1 

-1 n^ n^ 

f 8 = 1 8 = 1 

As in Section 6, the information that for each i and x the (j)i{x) are gaussian random 
variables with covariance given by Fq. (6.2), permits the averages over (j) in the 
definition of E5 and Eq to be evaluated analytically. It follows that E3 is given by 

E3 = Ee-—. (9.5) 



We then have 



E, = Es-2^E, + {^Y<[P-<P>,r>,. (9.6) 

6 6 
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Results obtained from 160 gauge field configurations with n^ of 140 are given in 
Table 5 and suggest that Li is a good approximation to tr log(A^^A^). 

A calculation of the effect of quark-antiquark vacuum polarization using a weak 
coupling perturbation expansion to leading order was reported recently in Ref. [9]. 
Staggered quarks are considered in Ref. [9] in place of our choice of Wilson quarks. 
The weak coupling expansion in Ref. [9] is expected to be reliable for sufficiently 
small gauge coupling and sufficiently large quark mass. For two flavors of quarks 
with degenerate mass ma in lattice units ranging from 0.05 to 1.00 and a parameter 
corresponding to /?' fixed at 5.68, the main effect of quark-antiquark vacuum polar- 
ization is found to be simply a coupling constant shift A/3. As ma ranges from 0.05 
to 1.00, the shifted /3 for full QCD runs from 5.34 to 5.63. 

For Wilson quarks, the mass is 

"^^ = ^-i- ^^-^^ 

Here kc is the critical hopping constant at which the pion mass becomes 0. With our 
choice of 0.1600 for k and with kc of 0.1694 [1] at (3' of 5.700, ma becomes 0.1734. To 
a first approximation, corresponding versions of QCD with Wilson quarks and with 
staggered quarks should have equal values of quark mass and (3. Thus the parameters 
of our trial calculation fall nearly within the range considered in the perturbative 
calculation of Ref. [9] and our results are qualitatively consistent with theirs. For 
two fiavors of staggered quarks with ma of 0.1734 and /?' of 5.68, the perturbative 
calculation predicts a A/3 of 0.226 in comparison to our prediction of 0.261 it 0.014 
for two fiavors of Wilson quarks. 

10 CONCLUSION 

The crucial question which we have not yet answered is how much time would 
be required to apply the algorithm we have described to QCD and determine, from 
Fq. (3.7) with n of 1, the error in the valence approximation to hadron propagators 
for more realistic choices of quark mass, lattice spacing and lattice volume than we 
chose in the test in the preceding section. If the algorithm can be run in reasonable 
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time with n of 1, it might be possible to use larger n and obtain smaller errors in 
hadron propagators. If the error found in this way for some small value of n is itself 
small for at least one quantity from which the physical value of the lattice spacing 
can be determined, it will be possible to calculate aj^ for full QCD. 

A perturbation theory estimate, which we will not discuss here, suggests that 
the optimal number of random (j) which our method requires will grow more slowly 
than a power of the inverse lattice spacing or the inverse quark mass. Similar es- 
timates suggest similar growth rates for the number of independent gauge configu- 
rations needed to evaluate the expectation values entering the determination of the 
coefficients a^ in the expansion in Eq. (3.1). The remaining question is how large an 
ensemble of gauge configurations may be required for small values of lattice spacing 
and quark mass to find the valence approximation error in hadron propagators using 
Eq. (3.7). If the difference Qn for some small value of n turns out to be quite small, 
as occurs for the parameter values in Section 9, or if Qn is sensitive only to low mo- 
mentum fluctuations of the gauge held, the calculation of propagator errors may be 
possible with reasonable ensembles sizes. We do not know at present whether one of 
these conditions might be realized for values of lattice spacing and quark mass which 
would permit an extrapolation to the physical limit of hadron masses. 
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A CONVERGENCE PROOE 

We now prove convergence for any bounded G of the sequence of approximate 
vacuum expectations < G >„ deflned by Eq. (3.6). To begin, we will show that 
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convergence holds if the remainders Qn in Eq. (3.5) have absolute value bounded by 
a constant independent of n and the gauge field. The required bound on \Qn\ we will 
then show follows from bounds on the coefficients of Qn expanded in the basis fi of 
Section 2, on the vectors fi multiplying these coefficients and on the number of terms 
occuring in this expansion for any fixed value of the dimension sum di. 
It is convenient to recast Eqs. (2.4) in the form 

<G> = I 

Y = C^ [ diiGexp{AS + S), 

Z = C^ f dfiexpiAS + S), (A.l) 

and Eqs. (3.6) in the form 

<G>n = — 

Yn = C jdiiGexpiAS-Qn + S), 

Zn = C fdiiexp{AS -Qn + S), (A.2) 

where S is the bounded effective action entering the inner product of Eq. (2.2), ( is 
defined in Eqs. (2.2), Qn is defined by Eq. (3.5) and AS is given by 

AS = logdet{M^M + R) + -P-S. (A. 3) 

The difference Y — Yn can be expressed in the form 

Y-Yn = C l'dfiGQn^^^^^^^^exp{AS + S). (A.4) 



The discussion of Sections 2 and 3 combined with the boundedness of S implies AS is 
bounded in absolute value. We will show below that Qn is bounded in absolute value 
by a constant independent of n and the gauge field. It follows that [1 — exp[ — Qn)]/Qn 
is bounded in absolute value by a constant independent of n and gauge field. Thus 
from Eq. (A.4) it follows that there is a constant c independent of n such that 

\Y - Yn\ < cC' Jdii \G\ \Qn\exp{S). (A.5) 
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By the Cauchy-Schwarz inequality we then have 

|i^-i;i<c||G'||||g„||, (A.6) 

where || • • • || is the norm defined by Eq. (2.1). The discussion of Section 3 imphes 
\\Qn\\ goes to as n — > oo. Thus |y — 1^1 goes to 0. The preceding argument with G 
chosen to be 1 imphes that \Z — Zn\ goes to 0. Thus by Eqs. (A.l) and (A. 2), < G > 
approaches < G >„. 

To show that Qn is bounded in absolute value by a constant independent of 
n and the gauge field, it is easier to work with a related remainder Q'^ found by 
expanding log det{M^M + K) using the basis fi of JF. In JF we have 

logdet{M^M + R) = ^6j„ (A.7) 

i 

h = {h,\ogdet{M^M + R)). (A.8) 

Define the remainder for truncations of this expansion Q',^ to be 

Q'n = Y.bJ.- (A.9) 

On 

Each hi in the basis for "H is a linear combination of fj in the basis for JF with a 
single fixed value of total dimension dj. In addition, distinct hi and hii are linear 
combinations of disjoint sets of fj. It follows that there is a way of choosing the 
details of the ordering of basis vectors fj which we specified in Section 2 which makes 
the remainder sequence Qn a subsequence of the remainder sequence of Q'^. It is 
therefore sufficient for us to establish the existence of a bound on \Qn\ independent 
of n and the gauge field. 

B EXPANSION COEEEICIENT BOUND 

As discussed in Section 3, since M is a finite matrix with matrix elements 
bounded uniformly over all gauge fields, the spectrum of M^ M + i? is bounded from 
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above by some constant A. The spectrum of M^ M -\- R is also bounded from below 
by R. We can therefore express log det{M^M + K) in the form 

logdet{M^M + R) = logdet{^^^-^ + ^) + dMlog{A), (B.l) 

M^M R 
= trlog[l-{l-^--j)] + dMlog{A). (B.2) 

where o?m is the dimension of the matrix M. The spectrum of the matrix 1 -r -j 

is nonnegative and bounded from above by 1 — ^ independent of the gauge fields. 
Thus the logarithm in Eq. (B.2) can be expanded as a power series, 

logdetiM^M + K) = -Y^-tr[il-^^-jy] + dMlogiA), (B.3) 

which converges as a result of the bound, independent of gauge field, 

M^M R R 

tr[{l-^-jy]<dM{l-jy. (B.4) 

The convergence bound Eq. (B.4) permits Eq. (B.3) to be substituted into 
Eq. (A. 8) to obtain an expansion for the coefficients 6i, 

3 

1 -^ M'^M R ■ 

b.3 = -jif^Mil-^-jY])- (B.5) 

Erom this expansion we obtain a bound on the bi. 

The matrix elements of M^ M are either independent of the gauge field or are 
linear functions of matrix elements of the product u[x^ y)u{y, z) for a pair of adjoining 

nearest neighbor links (x, y) and (y, z). Thus (1 -^ -jY includes products of at 

most j matrix elements of u{x^ y)u{y, z). Eor the trace tr[(l -^ ^)-'] written as 

a linear combination of the irreduceable representations fi discussed in Section 2, we 
will show that the dimension sums di which occur are bounded by 

J, <£ + i(j + 3)3, (B.6) 
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where £ is the number hnks in the lattice. 

Standard results on products of SU{3) representations imply that for a single 
link the largest dimension d of an irreducable representation which can occur in the 
product of pi copies of matrix elements of u[x^y) and p2 copies of matrix elements 
of the conjugate matrix u(y, x) corresponds to a Young tableau with pi columns 
containing one box and p2 columns containing two boxes. The dimension d is then 
given by 

d=^ipi + l)ip2 + l)ipi+P2 + 2). (B.7) 

Now consider the dimension sums di which occur if a product of j matrix elements 
of u[x^y)u[y^ z) for any collection of j pairs of adjoining links (x,y) and (y^z) is 
expressed as a sum of the products of irreduceable representations fi. Since d in 
Eq. (B.7) rises faster than linearly in pi and p2, the dimension sum di will be great- 
est if all j matrices u[x^y)u[y^ z) are taken on a single pair of adjoining links. By 
Eq. (B.7) this will be maximized if pi matrix elements of u[x^y)u[y^ z) are chosen 
and p2 elements of the conjugate matrix u[z^ y)u{y, x) are chosen with pi given by the 
largest integer less than or equal to j/2 and p2 is given by j — pi. The contribution to 
di coming from these two links is 2d for d given by Eq. (B.7). This sum is bounded 
by the second term in Eq. (B.6). The remaining lattice links all carry the trivial 
1-dimensional representation giving a sum bounded by the first term in Eq. (B.6). 

The /i, however, have been ordered with increasing values of o?i, and the fi 
come from a Gram-Schmidt process applied to the fi. Thus fii is orthogonal to any 
fi with di less than dii . Therefore hi^ in Eq. (B.5) vanishes unless j obeys Eq. (B.6). 
The Cauchy-Schwarz inequality applied to Eq. (B.5), combined with Eq. (B.4) and 
(B.6) then gives the bound 

\h\<c{l-^f''''\ (B.8) 

for a pair of constants c, c' independent of i. Eor the infinite set of i large enough 
that di is much larger than £, the | power in Eq. (B.8) follows from the power 3 in 
Eq. (B.6). By choosing c sufficiently large, Eq. (B.8) can then be made to hold also 
for the remaining finite set of smaller i and di. 
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C BASIS VECTOR BOUND 

We next derive a bound on the fi independent of the gauge field. This bound 
is obtained by examining the Gram-Schmidt process by which the sequence of fi is 
constructed from the sequence of fi. 

The vector fi can be found by normahzing the vector pi 

h = Z^, (c.i) 

Pi = f^ + Y^^'^jfj^ (*^-2) 

with coefficients aij, defined only for j less than z, chosen to minimize the norm 

Ib.ll^ = (/. + Yl "u/j. /«■ + Yl "u'/jO- (C-3) 

By the triangle inequality and Eq. (C.2), we have for any gauge field and any i 

The fj are normalized with respect to pure Haar measure or, equivalently, using 
the inner produce of Eq. (2.2) with S replace by 0. Eor any irreducable representation 
of SU{3) by unitary matrices Vab^ of dimension d X d^ we have for an integral over 
Haar measure i/ on a single copy of SU{3) 

v*^Vabdu = -. (C.5) 

Thus the matrix ydvab is normalized to 1 and has matrix elements bounded by yd. 
Applying this argument to the copy of SU{'i) on each lattice link and then taking a 
maximum of the resulting product of dimensions we obtain the bound, independent 
of gauge field, 

l/.l<(f)^- (C.6) 
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Eqs. (C.2) and (C.6), combined with the Cauchy-Schwarz inequahty and the infor- 
mation that the sequence of dj is nondecreasing then imply for any gauge field 

b.l<(i + EK-r)'(7)'- (C-7) 

Now return to the squared norm (pi^pi). Since S in Eq. (2.2) defining the norm 
is bounded from above and below, and since the fj are orthonormal with respect to 
Haar measure, we have for a constant c independent of i 

iP^,Pi) = C^ j d^\f^ + ^a,Jj\^exp{S), (C.8) 

3<i 

> cJdi,\f, + J2a,,f,\' (C.9) 

= c(l + ^|a.,f). (CIO) 

Eqs. (C.l), (C.7) and (C.8) then yield for some constant c' independent of i and the 
gauge field 

l/.l<c'(|)i (c.ii) 

D REMAINDER BOUND 

To complete the bound on the remainders Q'^ we need a bound on the number 
n[d) of fi entering Eq. (A. 9) for any fixed value d of the dimension sum di. The 
number n[d) of such fi is certainly bounded by 

n{d) < [d^m{d)Y, (D.l) 

where m[d) is the number of distinct irreducable representation of dimension d or 
less for a single copy of SU{3) and d^ is clearly an upper bound on the number of 
matrix elements in each such matrix. Each irreduceable representation of dimension 
d or less corresponds to a Young tableau specified by pi and p2 as before but with 
Eq. (B.7) now replaced by an inequality 

d>-ipi + l)ip2 + l)ipi+P2 + 2). (D.2) 
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Any pi or p2 fulfilling Eq. (D.2) is less than v 2o?. Thus the total number of such pairs 
is less than 2d. Eq. (D.f) becomes 

n{d) < {2dy. (D.3) 

Combining Eq. (A. 9) for the remainder Qn with the bounds of Eqs. (B.8), 
(C.ii) and (D.3) we obtain a bound of the form 



dydn 



A' 



\Q'n\< E^"^^(1-t)^^'''^- (D-4) 



The sum in Eq. (D.4) is over all integers d greater than or equal to the dimension sum 
dn for /„. The parameters i?, A, c and c" are independent of d and the gauge field. 
Since dn is a positive nondecreasing function of n, the sum in Eq. (D.4) is less than 
the sum with o?„ replaced by 0, which in turn is convergent and gives the required 
bound on Q'^ independent of both n and the gauge field. 

E CHEBYSHEV EXPANSION 

We now derive the coefficients rii needed for the Chebyshev expansion of log y 
Eq. (6.7). Erom standard results on Chebyshev polynomials [13], it follows that 

-[1 + pT:^ii-r^)] = E ckm-—^), (E.l) 

y 1 e ^^Q 1 e 

where for 1 < k < n 

^ (1 + cosh x) sinh[(n + 1 - k)x] ,^ ^ 

sinh X cosh[(n + l)x] 



and in addition 



(1 + cosh x) sinh[(n + l)x] 
sinh X cosh[(n + l)x] 
-1 



cosh[(n + l)x] 
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(E.3) 



with X defined by 



coslix = . (E-4) 

1 — e 



Eq. (E.l) can be integrated from y to f using, for k greater tlian 1, tlie relation 



ni^y^ - ^ - tinr- (E-5) 



along with 



T^[x)dx 



t;{x) t*{x) 



8 8 



■> 



We obtain 



T*{x)dx = IliEl + Io{^_ (E.6) 



logy = 5: 6.r;(i-^) + eulogy (E.7) 

k=o -*- ^ 



(E.8) 



where for 1 < k < n 






hk 


(1- 


e)(l + cosh x) cosh[(n + 1 — h)x 




A;cosh[(n + l)x 


and in addition 








&n+i = 


(l-e)(l + coshx) 




2(n + l)cosh[(n + l)x]' 




&o = 


n+l 

k=i 



(E.9) 

The bound Eq. (6.8) on 5 follows from the integral of Eq. (E.l) combined with 
Eq. (E.2) for p and the bound 

mx)\ < 1, (E.io) 

for < x < 1. 
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O kg < \J Xmax I ^mtn > work/sweep work 

M^M 58.1 1.0 58.1 

N^N 0.091 22.1 2.0 44.2 

Table 1: Comparison of work required to evaluate \og{0)(f) for different choices of O 



rich El E2 

10 390.22 6870 

20 412.25 2733 

50 412.74 2843 

100 412.79 2848 

200 412.79 2848 



Table 2: Expectation values found from 16 gauge configurations each with 20 random 
(/), for various choices of the number of iterations of the Chebyshev algorithm used in 
the calculation of log{N^ N)(f). 
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n4, 


El 


E2 


A/3 


10 


416.8 ±2.1 


2340 ± 508 


0.247 ± 0.049 


20 


416.0 ±1.7 


2343 ± 387 


0.248 ±0.029 


30 


415.8 ±1.5 


2483 ± 367 


0.262 ±0.028 


40 


415.4 ±1.4 


2531 ± 343 


0.267 ±0.026 


50 


414.7 ±1.3 


2517 ± 344 


0.266 ±0.024 


60 


414.8 ±1.3 


2499 ± 318 


0.264 ±0.023 


70 


414.9 ±1.2 


2484 ± 309 


0.262 ±0.020 


80 


414.7 ±1.2 


2470 ± 294 


0.261 ±0.018 


90 


414.7 ±1.1 


2420 ± 290 


0.256 ±0.017 


100 


414.8 ±1.1 


2483 ± 307 


0.262 ±0.017 


110 


414.8 ±1.1 


2523 ± 309 


0.267 ±0.015 


120 


414.5 ±1.0 


2501 ± 296 


0.264 ±0.014 


130 


414.2 ±1.0 


2489 ± 290 


0.263 ±0.014 


140 


414.5 ±1.0 


2472 ± 288 


0.261 ±0.014 



Table 3: Expectation values and A/3 found from 160 gauge configurations for various 
choices of the number of random (f) used in the evaluation of traces. 
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loop 





valence 


valence + error 


full QCD 


n = 


n = 1 


n = 1, m = 1 


l3 = 5.439 


l3 = 5.439 


fi' = 5.700 


l3' = 5.700 


K = .16 



1 X 1 


0.4827 


;o5) 


0.5501 


(04) 


0.5501 


(29) 


0.5527 


:ii) 


2 X 1 


0.2438 


;o6) 


0.3261 


(06) 


0.3325 


(44) 


0.3341 


16) 


6i 


0.2747 


;o6) 


0.3618 


(05) 


0.3655 


(38) 


0.3678 


16) 


62 


0.2212 


;o6) 


0.3143 


(06) 


0.3208 


(42) 


0.3230 


:i7) 


3 X 1 


0.1243 


;o5) 


0.1966 


(06) 


0.2036 


(46) 


0.2064 


:i5) 


2x2 


0.0678 


;o5) 


0.1341 


(06) 


0.1479 


(46) 


0.1471 


16) 


4 X 1 


0.0635 


;o4) 


0.1188 


(05) 


0.1249 


(43) 


0.1283 


:i4) 


5 X 1 


0.0326 


:o3) 


0.0720 


(04) 


0.0774 


(35) 


0.0799 


10) 


3x2 


0.0195 


:o3) 


0.0592 


(04) 


0.0713 


(32) 


0.0706 


13) 


4x2 


0.0057 


:o2) 


0.0264 


(03) 


0.0346 


(24) 


0.0350 


;o9) 


3x3 


0.0038 


:o3) 


0.0206 


(03) 


0.0307 


(26) 


0.0288 


;o9) 



Table 4: Vacuum expectation of various Wilson loops found by four different methods. 



< [tr[log{N^N)]- < tr[log{N^N)] >i 
< [trlog{N^N)-Li]^ >i 



>i 



124 ±16 
16.0 ±6.3 



Table 5: Expectation values measuring how good an approximation Li is to 
tr log(A^^A^). The calculation uses 160 gauge configurations each with 140 random (/>. 
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A 

Q. 
O 

_o 

V 



1.0 



0.1 - 



0.01 



0.001 



4x15x13x24x23x3 



^ 



° D 



A ^ 



1x1 2x1 6i 62 3x12x2 



Figure 1: Vacuum expectations of various Wilson loops found by four different meth- 
ods. Boxes represent n = with (3 = 5.439, triangles are the valence approximation 
n = 1 with /?' = 5.700, circles are the valence approximation plus first error estimate, 
equivalent to n = 1^ m = 1 with /?' = 5.700, plusses are full QCD with (3 = 5.439 and 
two quark fiavors with k = O.fGOO. 
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